close all

dipole_data = load('C:\Users\PiascikSammy\Dropbox\QDM2_data\Sammy\MFZ-N1\MFZN1-A\FOV7\Gelson_code_boxes\results.mat', 'fits').fits;
bz_map = load('C:\Users\PiascikSammy\Dropbox\QDM2_data\Sammy\MFZ-N1\MFZN1-A\FOV7\2x2Binned\Bz_uc0.mat', 'Bz').Bz;

dipole_data.log_moments = log10(dipole_data.moments);
gelson_coords = readmatrix('C:\Users\PiascikSammy\Dropbox\QDM2_data\Sammy\MFZ-N1\MFZN1-A\FOV7\FOV7_sources.csv');
box_widths = gelson_coords(:, 2) - gelson_coords(:, 1);
box_heights = gelson_coords(:, 4) - gelson_coords(:, 3);
gelson_coords = [gelson_coords(:, 1), gelson_coords(:, 3), box_widths, box_heights];

right_direction = abs(dipole_data.decs-180)<90;
near_surface = dipole_data.heights > -1.5e-5;

magnetite_density = (5.17/1000)*(100)^3; % In units of kg/m^3
magnetite_saturation = 83; % In units of Am^2/kg

% Naive assumption of cubic magnetite grains
dipole_data.lengths = (dipole_data.moments/magnetite_saturation/magnetite_density).^(1/3);
big_enough = dipole_data.lengths > 2e-7;

QDM_figure(bz_map);
hold on
targets = gelson_coords(right_direction & big_enough & near_surface, :);
for i = 1:size(targets, 1)
    rectangle(Position=targets(i, :));
end

figure
histogram(dipole_data.heights(right_direction & big_enough))
